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ABSTRACT 

For the very first time we present 3D simulations of planets em¬ 
bedded in stellar irradiated discs. It is well known that thermal 
effects could reverse the direction of planetary migration from 
inwards to outwards, potentially saving planets in the inner, op¬ 
tically thick parts of the protoplanetary disc. When considering 
stellar irradiation in addition to viscous friction as a source of 
heating, the outer disc changes from a shadowed to a flared struc¬ 
ture. Using a suited analytical formula it has been shown that in 
the flared part of the disc the migration is inwards; planets can 
migrate outwards only in shadowed regions of the disc, because 
the radial gradient of entropy is stronger there. In order to con¬ 
firm this result numerically, we have computed the total torque 
acting on planets held on fixed orbits embedded in stellar irra¬ 
diated 3D discs using the hydrodynamical code FARGOCA. We 
find qualitatively good agreement between the total torque ob¬ 
tained with numerical simulations and the one predicted by the 
analytical formula. For large masses (> 20Mq) we find quanti¬ 
tative agreement, and we obtain outwards migration regions for 
planets up to 6 OM 0 in the early stages of accretional discs. We 
find nevertheless that the agreement with the analytic formula 
is quite fortuitous because the formula underestimates the size 
of the horseshoe region; this error is compensated by imperfect 
estimates of other terms, most likely the cooling rate and the 
saturation. 


1 INTRODUCTION 


In recent years it has been shown that low mass 


discs with non-isothermal effects (Paardekooper 

& Mellema 2006; Baruteau & Masset|2008j Kley 

& (Jrida 2008 Kley et al. 2009, 

Masset and Ca- 

soli 2009 

Masset & Casoli 2010; Paardekooper et 

al. 2010 

2011 Lega et al. 2014 

). Precisely, the 


migration in the inner part of a radiative disc can 
be directed outwards, while it remains directed 


inwards in the outer disc (Bitsch and Kley 2011). 


This establishes the existence of a critical radius 
where migration vanishes, towards which plane¬ 
tary cores migrate from both the inner and the 
outer part of the disc. Therefore, the zero mi¬ 


gration location acts as a planet trap at which 
proto-planets can accumulate in resonances, col¬ 


lide and eventually form bigger objects. (Lyra et 


aL]| 2010| ICossou et al.||2014| |Coleman fe Nelson 

2014| >. 


In all these works the heating is provided 
by viscous friction and the cooling by radia¬ 
tive diffusion (in 3D discs) or by a local cool¬ 
ing rate (in 2D discs). More recently, 
al. (2013]), have shown the importance 
irradiation on the disc structure and the conse¬ 
quences for planetary migration. The main re¬ 
sult is that when considering stellar irradiation 
the outer disc changes from a shadowed to a 
flared disc. Nonetheless, opacity transitions (e.g. 
the iceline) create bumps in the aspect ratio and 


Bitsch et 
of stellar 
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hence shadowed regions. In the flared part of the 
disc the migration turns out to be inwards and 
planets can migrate outwards only in shadowed 
regions of the disc, where the radial gradient of 
entropy is much steeper than in the flared part 
of the disc. These results were obtained in equi¬ 
librium discs with uniform viscosity featuring 
a zero radial mass flux ( Bitsch et al.p 013l and 
extended to the case of accretion discs with an 


GOCA (Lega et al. 2014) that includes a two- 


alpha prescription for the viscosity (Shakura & 
Sunyaev]|1973 1 in Bitsch et al. (2014 1 . This sec¬ 


ond case is of particular interest, since a given 
mass-flow (M) rate corresponds to a specific disc 
age (Hartmann et al. 19981 so that, investigat¬ 


ing the disc structure for different values of M, 
is equivalent to studying the disc structure as a 
function of the disc evolution (Bitsch et al. |2014| 
2015). 


In |Bitsch et al.| ( f2013| |2014[|2015| ) the planet 
migrations maps have been obtained applying 


the torque formula of Paardekooper et al. (2011) 


using the disc properties obtained in the sim¬ 
ulations. The formula, although very complex, 
is based on the fact that the torque exerted 
by the protoplanetary disc onto the planet has 
two main contributions: i) the so-called Lindblad 
torque due to the spiral arms launched by the 
planet in the disc, which is not affected by the 
equation of state 0 and ii) the co-orbital coro¬ 
tation torque caused by material librating in the 
horseshoe region. When considering radiative ef¬ 
fects the corotation torque contribution can be 
positive and possibly dominate over the negative 
Lindblad torque, leading to outwards migration. 
The Paardekooper et al. formula was calibrated 
with a 2D hydrodynamical model for low mass 
planets (5M0). 

With a specific disc setting not accounting 


for stellar irradiation, Kley et al. (2009) found 


positive total torque for planetary masses in 
the range [5,30] Earth masses. The influence of 
the disc’s mass on the migration was studied in 
Bitsch and Kley] (2011) and its application to 
Earth sized planets was addressed by |Lega et] 
al. (2014) who found the contribution of a new 


torque not accounted for by the formula. This 
new negative torque is due to the formation of an 
asymmetric cold and dense finger of gas driven 
by circulation and libration streamlines. 

The aim of this paper is to investigate nu¬ 
merically the total torque acting on planets kept 
on fixed orbits in 3D stellar irradiated discs 
and provide quantitative comparisons with the 


Paardekooper et al. (2011) formula. We use 


temperature solver for radiative transfer in the 
flux-limited approximation. 

In principle stellar heating only changes the 
structure of the disc and does not act on the 
mechanism responsible for outwards migration 
directly. Nevertheless the analytical formula has 
been calibrated on 2D discs so that a quanti¬ 
tative test on the validity of that formula for 
realistic 3D discs is needed. Moreover, a con¬ 


stant mass flow M through the disc (Bitsch et al. 


20141 might have an influence on the torque act¬ 


ing on planets and, to our knowledge, this case 
has never been tested before in 3D discs with a 
non-isothermal EOS. 

The paper is organized as follows: the physi¬ 
cal modelling is presented in Section 2, in Section 
3 we describe the migration maps obtained from 
the analytic formula applied to our 2D unper¬ 
turbed discs [^] In Sections 4 and 5 we provide 
results of 3D simulations done respectively on a 
constant viscosity stellar irradiated equilibrium 
disc and on an a-viscosity accretion disc. In sec¬ 
tion 6 we go beyond the simple comparisons of 
the net torques predicted analytically and mea¬ 
sured numerically, focusing on a key ingredient 
of the analytic formula: the size of the corota¬ 
tion zone, which governs the corotation torque 
and the torque saturation. The conclusions are 
provided in Section 7. 


2 PHYSICAL MODELLING 

The protoplanetary disc is treated as a three di¬ 
mensional non self-gravitating gas whose motion 
is described by the Navier-Stokes equations. We 
use spherical coordinates (r, 9, ip) where r is the 
radial distance from the star, i.e from the ori¬ 
gin, 9 the polar angle measured from the 2 -axis 
(the colatitude) and ip is the azimuthal coordi¬ 
nate starting from the a:-axis. The midplane of 
the disc is at the equator 9 = f. We work in 
a coordinate system which rotates with angular 
velocity: 


G(AL* -(- nip) 


GM* 


the explicit/implicit hydrodynamical code FAR- 


where M* is the mass of the central star, G the 
gravitational constant and a p is the semi-major 
axis of a planet of mass m p , assumed to be on a 
circular orbit. The gravitational influence of the 
planet on the disc is modelled as in | Kley et al.| 
|2009) using a cubic-potential of the form: 


Actually, it scales with yr 0 , with Fq given in Eq[| 


2 discs in thermal equilibrium in the (r — z) plane 
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d > e 
d < e 


(1) 


where d is the distance from the disc element to 
the planet and e is the softening length. Writing 
x = d/e the function / is given by: 

f{x) = x 4 — 2x 3 + 2x . 

We have considered e = 0.6 Rh in our sim¬ 
ulations sets SED and SAD and t = 0.5Rh for 
simulation set RED, with Rh the Hill radius of 
the planet: 


Rh = a f 


3 M t 


We will discuss in Section 6 the dependence of 
the total torque on the softening length. 

In the Paardekooper formula, there is also 
a softening length parameter f3 and it has been 
fixed to /3 = 0.35 H (H being the disc’s local 
scale height) from previous comparison between 
the formula and 3D simulations (Bitsch and Kley 


2011 ). 


The hydrodynamical equations solved in the 


code are described in Lega et al. (20141, the two 


temperatures approach for the stellar irradiation 


was described in detail in Bitsch et al. (2013) and 


the opacity prescription in Bitsch et al. (20141. 


The flux received from the star at a radial dis¬ 
tance r is: 


F* = RfaTf/r 2 


( 2 ) 


were the radius of the star is set to R* = 1.5Rq 
and the temperature T* = 4370K (appropriate 
for a Solar-type protostar). We consider the disc 
settings of (Bitsch et al. 2013) namely a disc ex¬ 
tending from r min ^ r ^ r max with r m in = 1AU 
and rmax = 50AU. In the vertical direction the 
disc extends from the midplane (9 ~ 90°) to 20° 
above the midplane, i.e. 9 ~ 70°). The initial 
surface density profile is E(r) = Eo (r/aj)~ b and 
aj = 5.2AU. 

We consider in the following different sets 
of simulations. Precisely, we consider a stellar ir¬ 
radiated equilibrium disc (set SED, hereafter), 
with uniform viscosity (viscosity coefficient v = 
10 _ 5 ajllp), Eo = 4.88 x 10 ~ 4 in code units 
(147g/cm 2 ) and b = 0.5. We also consider a 
stellar irradiated accretion discs (set SAD, here¬ 
after) with alpha viscosity (a = 0.0054), con¬ 
stant M = 4 ■ 10 -8 MQ/yr and initial value of 
b = 15/14 and of Eq = 430g/cm 2 . The ini¬ 
tial value of b corresponds to the radial sur¬ 
face density profile in flared disc [^] with con¬ 
stant M at all orbital distances. A disc with 
M itm 4 ■ lO~ 8 M 0 /yr corresponds to the early 


3 with flaring index 2/7 
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evolution stages of accretion discs and is proba¬ 
bly younger than 1 Myr ( Hartmann et al.|1998 l. 
The disc structure of the SAD case is discussed 


Bitsch et al. (2015 I. 


For comparison, we will also consider the 
case of a radiative equilibrium disc with uniform 
viscosity (v = 10 - 5 a 2 f2 p ) and b = 0.5, were 
the only source of heating is viscous heating (set 
RED hereafter). This setup of the RED disc 
has been previously studied in Kley et al. ( |2009| ; 
|Bitsch and Kley] ( |2011[ ) and in |Lega et al'| ‘ |2014[ ) 
for planetary masses up to 3 OM 0 . 

Before placing the planet in a 3D disc, we 
bring the disc to radiative equilibrium. We first 
model each disc in 2D, with coordinates (r, 9). 
For the accretion disc the (r, 9) disc evolution is 
explained in detail in Bitsch et al. (2014 |2015| ). 
Once the 2D equilibrium is achieved all the gas 
fields are expanded to 3D. The aspect ratio of 
the 2D equilibrium for the SED, the SAD with 
M = 41O _8 M0/yr and the RED are reported 
in Fig[l] We remark that for the SED and RED 
with uniform viscosity, the settling of the aspect 
ratio to an equilibrium value does not change the 
surface radial density profile of the disc. Instead, 
in the SAD disc, in which the viscosity follows 
the alpha-prescription and thus depends on the 
disc’s local scale height H, the resulting surface 
density profile at equilibrium is somewhat differ¬ 
ent from the initial one (see Bitsch et al. (j2014|). 

In the two cases where the disc is heated by 
the star we observe the typical profile of a flared 
disc, while we have a shadowed disc in the RED 
case. The bump in H/r around 3—4 AU is caused 
by a transition in opacity that changes the local 
cooling rate of the disc and hence changes the 
disc’s temperature (Bitsch et al. 2014 1 . 

The resolution of our computational grid is 
chosen in order to have in the radial direction 
approximatively n grid-cells in the horseshoe re¬ 
gion. The half-width of the planet’s horseshoe 
region is given, in the isothermal disc approxi¬ 
mation (Masset et al.|20'06|, by: 


Xhs 2 D = 1.16a p * f 


( 3 ) 


were q = m p /M * and h is the disc aspect ra¬ 
tio at a p . We have checked that this formula 
can still be used as a guideline for the choice of 
the resolution for non isothermal 3D discs with 
and without stellar irradiation. However, Eqj3] 
is based on the equivalence of the linear corota¬ 
tion torque and the horseshoe drag in 2D discs. 


In Masset et al. (2006) it has been shown that 


nonlinearities appear on the flow for mass ratios 
q ^ h 3 . For these planetary masses the horse¬ 
shoe region is larger than the linear prediction 
resulting in a boost of the corotation torque. It is 
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Figure 1. Aspect ratio of the 2D (r, 6) equilib¬ 
rium of simulation set SED (green) and SAD with 
M = 410“ 8 (blue). For comparison the red curve 
provides the aspect ratio of the radiative equilibrium 
disc (RED) studied in|Kley et al. | (|2009|l ; |Lega et al.| 
( |2014| ) . 


therefore important to check whether also in 3D 
hydrodynamical simulations the horseshoe width 
increases for q ^ h 3 and what is the correspond¬ 
ing impact in the corotation torque. Section 6 is 
devoted to this topic. 

The values of the masses and resolutions 
(N r , Nq, Nip) for the SED disc are shown in Ta¬ 
ble [I] In Lega et al. ( 20141 we have shown that 
a resolution of 4 grid cells in the half width of 
the horseshoe region matched the requirement of 
having resolution independent results within rea¬ 
sonable CPU times for masses up to 2 OM 0 . Here 
we use 4 grid cells in the half width of the horse¬ 
shoe region also for planets with masses larger 
than 20M ffi ; We have tested that results are sta¬ 
ble when increasing the resolution to 5 grid cells 
in some test cases. 

Concerning the calculation of the gravita¬ 
tional torque acting on the planet we recall that 
it is common to exclude the inner part of the Hill 
sphere of the planet. This is obtained by apply¬ 
ing a tapering function (named Hill cut in the 


following) which in our case reads (Kley et al. 


2009): 


Md) = 


exp 


d/R h — 6 

6/10 


+ 1 


(4) 


The value of fb is zero at the planet location 
and increases to 1 at distances d from the planet 
larger than the Hill radius Rh- The parameter 
b denotes the distance from the planet, in unit 
of Hill radius, at which fb = 1/2. Here we use 
b = 0.8 ( |Crida et al.| [2009). 


This procedure allows to exclude the part of 
the disc that is gravitationally bound to planets 
that form a circum-planetary disc. However, this 
prescription is not justified for small mass plan¬ 
ets that do not have a circum-planetary disc. For 
small planets we mean planets having a Bondi 
radius (Rb) smaller than their Hill radius. From 
the definition of the Bondi and Hill radius one 
obtains Rb < Rh for q < h 3 /%/3 |^j In the 
following we show the total torque computed 
with Hill cut, and, for planetary masses satis¬ 
fying q<h 3 /\/ 3, we also show the total torque 
computed without Hill cut. 


3 MIGRATION MAPS 


The change in the disc structure due to stel¬ 
lar irradiation has important consequences for 


the migration of embedded bodies (Bitsch et 


al. 20131. One can estimate the torque act¬ 


ing on planets using the formula provided by 


Paardekooper et al. (2011). The migration maps 


obtained from the formula are shown on Fig[2]for 
respectively set SED (top panel), set SAD with 
M = 4 lO _ 8 M 0 /yr (middle panel) and set RED 
(bottom panel). The values of the total torque 
Ttot are are normalized with respect to 


To = (q/h) 2 X p a 4 p n 2 p 


(5) 


where E p is the disc’s surface density at the 
planet location a p . 

We do not enter here in the details of the 
torque formula provided by |Paardekooper et al.[ 
(201 lj) we just recall that outwards migration is 
a delicate process depending on the dynamical 
properties of the corotation region as well as on 
the viscosity of the disc and of its cooling proper¬ 


ties. Different timescales are at play (see Bitsch 
& Kley] (2010}), in order to detect a positive coro¬ 
tation torque contribution which possibly domi¬ 
nates over the negative Lindblad torque. 

We remark that the formula has been ob¬ 
tained for planets that do not perturb the disc 
significantly (q < h 3 , i.e. linear regime) while the 
maps shown in Fig [2] consider also intermediate 
mass planets (30 — 70 M^). The torque satura¬ 
tion, which determines the upper limit in planet 
mass for outwards migration, depends on the vis¬ 
cous timescale defined by: 


r v = (xhsZD'Y 1/4 ) 2 /^ 


( 6 ) 


were 7 is the adiabatic index, 7 = 1.4 in our 
simulations. The width of the horseshoe region 


4 We notice that the parameter that determ ines the 
flow linearity in the planet vicinity is q < h 3 ( |Masset| 
|et al.|2006[ ). 
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5 


is Xhs 2 D / 7 1 in order to take adiabatic effects 
into account (Paardekooper & Papalo izou[200 9). 
However, for intermediate mass planets, the ac¬ 
tual width of the horseshoe region can be dif¬ 
ferent from Xhs 2 D of Eqj3] (Masset et al. 2006) 
and this difference can have an impact on torque 
saturation. 

I 11 the migration map of the SAD case 
(Fig|jtop) we see that, using Xhs 2 D / 7 1 / 4 , torque 
saturation is supposed to occur for quite large 
planetary masses and outwards migration is ex¬ 
pected for planets up to 70 Mq (this upper limit 
can moderately change by changing the smooth¬ 
ing length /?). 

It is interesting to check this result with 3D 
numerical simulations, because outward migra¬ 
tion of intermediate mass planets has important 
impacts for formation models of giant planets. In 
fact giant planet precursors could be prevented 
from migrating into the inner disk before they 
reach a mass that allows gap opening and there¬ 
fore slower inward migration in the type-II mi¬ 
gration regime. 

I 11 order to test the validity of the formula we 
have considered planets of different masses held 
on fixed orbits in each of the three considered 
discs and we have computed the evolution of the 
torque with time until we obtain a stationary 
state. We choose the distance of the planets from 
the star within the region of expected outwards 
migration from Figj2] 


4 EQUILIBRIUM DISCS 


In the SED we have embedded planets of differ¬ 
ent masses held on a fixed orbit at r = 4AU, for 
which we expect outwards migration from the 
torque formula (Fig[2] middle panel). Let us no¬ 
tice that, in the inner part of the disc, the as¬ 
pect ratio of the SED is very similar to that of 
the RED (Fig|T]y ; since the density gradient is 
the same for both discs the computation of the 
total torque should give results similar to those 


obtained in Kley et al. (20091 and Lega et al. 
(2014). 


When considering planets at 4AU we ob¬ 
serve (Fig[3| that for masses larger than 30 Earth 
masses the torque is slightly larger than the 
one provided by the formula and the transition 
to negative torque occurs at 45 M 0 instead of 
4 OM 0 as expected from the formula (Figj2j mid¬ 
dle panel). However, the overall picture is quite 
in good agreement with the results provided by 
the analytical formula. For comparison we have 
extend our previous study of the RED from |Lega| 
et al. (20141 for planets at 5.2AU towards plan¬ 


ets with masses up to 60 Earth masses. We can 
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Figure 2. Migration maps for simulation set SAD 
with M = 410 -8 MQ/yr (top), SED (middle). For 
comparison, the migration map of the radiative equi¬ 
librium disc (RED) is also plotted (bottom). The 
white x symbols indicate the position and the masses 
of the planets for which we compute the total torque 
with 3D simulations in this paper. 
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mass (Mq) (N r ,Ng, N v ) n cells in Xh s 2 D 


8 

(1326,66,908) 

4 

10 

(1186,66,812) 

4 

20 

(1086,82,740) 

5 

25 

(1010,66,702) 

4 

25 

(1238,80,876) 

5 

30 

(911,66,620) 

4 

40 

(770,66,542) 

4 

50 

(670, 60, 468) 

4 

60 

(612,60,430) 

4 

60 

(764,66,540) 

5 


Table 1. Simulations parameters for set SED. Embedded planets are at dis¬ 
tance 4AU from the star. 


appreciate in Fig|4] the agreement between the 
results of numerical simulations and the values 
given by the analytical formula. 

In both the RED and SED, we provide in the 
following subsections a quantitative comparison 
between the values of the torque resulting from 
numerical simulations and the values given by 
the formula. 


4.1 Planetary masses smaller than 

1OM 0 


In the case of the SED disc, the minimum 
planetary mass to have outwards migration is 
somewhat larger than that predicted in the 


Paardekooper et al. (2011) formula ( about 8M® 


against 4M®). This minimum mass is very im¬ 
portant for giant planet formation models. The 
larger transition mass that we find can be ex¬ 
plained with the cold finger effect studied in |Lega| 
[et al.| ( |2014| ). The unperturbed temperature at 
the planet location is of 112K, and, according to 
the criterium given in Lega et al. (2014), a nega¬ 
tive contribution to the total torque (cold finger) 
exists when the Bondi radius Rb is smaller that 
the Hill radius Rh- We have suggested: 


Rb = K/7.1M $ ) 2 ^ 3 

Rh (a/5.2AU)(T/75K) ’ l ' 

Applying Eq[T] for planets at 4AU where T = 
112 K, we obtain = 1 for a 1OM 0 and 
=0.6 for a 5M 0 planet. Thus, the expected 
transition between outwards and inwards migra¬ 
tion occurs in the interval [5,1O]M 0 in very good 
agreement with the value of 8M® found in the 
simulations. In the RED disc, the minimum 
planetary mass to have outwards migration is 
about 4M 0 again st 2M 0 , due to the cold fin¬ 
ger effect found in |Lega et al.| (|2014|) . 


4.2 Planetary masses in the interval 

10 - 3OM 0 

In Fig[4] we notice that the maximum value of 
the computed torque as a function of planetary 
mass is about a factor 2 smaller than the max¬ 
imum provided by the formula. The planetary 
mass for which we measure the largest positive 
torque value is respectively of 2OM 0 for the RED 
case and of 3OM 0 for the SED case. In both cases 
the formula provides the largest positive torque 
for a planet of 1OM 0 , the onset of saturation 
occurs at planetary masses smaller than what 
found in 3D simulations. We will discuss this 
point in Section 6. These quantitative discrep¬ 
ancies have already been pointed out by |Kley ~et| 
al. (20091 caused by differences between 2D and 


3D effects. 


4.3 Planetary masses larger than 3OM 0 

The analytical formula was derived for small 
mass planets (5M 0 ), i.e. for planets that only 
slightly perturb the disc. However, larger mass 
planets start opening a partial gap in the disc, 
meaning they significantly perturb the disc. 
These perturbations make the application of the 
torque formula dubious for masses larger than 
3OM 0 . However, the density in the corotation re¬ 
gion is perturbed by less than 30% up to 60M®; 
i.e we are still far from planetary masses which 
really open a gap in the disc. 

We remark that the formula and simula¬ 
tions results are in very good agreement for 
masses larger than 3OM 0 . The torques in the for¬ 
mula are calculated from the gradients of surface 
density, temperature and entropy in the unper¬ 
turbed disc. In the formula, the corotation torque 
saturates as the planet mass increases, making 
the total torque transit to negative. In our simu- 
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Figure 3. Stellar irradiated equilibrium disc (SED): 
comparison with Paardekooper et al. ( 2011 ) formula 
for planets on fixed orbits at 4AU. For planets of 
10 and 8 M 0 we also show the value of the torque 
computed without the Hill cut (see text). 



m p (Earth masses) 


Figure 4. Same as Fig.3 for the radiative disc (RED) 
studied in Lega et al.| ( |2014[ ) whose aspect ratio is re¬ 
ported in Fig[l]and the corresponding migration map 
in Fig[2] bottom panel. For planet of mass smaller 
than 1 OM 0 we also show the value of the torque 
computed without the Hill cut (see text). 


lations, the total torque is obtained directly from 
the density structure of the disc, meaning that 
on top of saturation effects, the partial opening 
of a gap is also taken into account, decreasing 
even more the corotation torque. That these two 
different approaches match in the total torque 
may be a coincidence, but we remind that the 
density in the corotation region is perturbed by 
less than 30% up to 60M®. 

We can conclude this section saying that 


when we take into account stellar irradiation in 
an equilibrium disc no additional effects on the 
total torque acting on planets in the explored 
mass range of [8, 60] Earth masses are observed. 
Moreover, the analytical formula Paardekooper] 


et al. (2011) can be used with good confidence 


when performing, for example, N-body simula¬ 
tions with a migration prescription. 


5 STELLAR IRRADIATED 
ACCRETION DISCS 


Protoplanetary discs accrete gas onto the central 
Star (Lynden-Bell & Pringle 19741. Thus equi¬ 
librium discs like the SED and RED cases in¬ 
vestigated above are not a good approximation 
of reality, given that they do not transport mass 
radially. A much better approximation are the 
discs where the mass flow M is independent of 
radius, called accretion discs. Their structure has 
been investigated in |Bitsch et al. (20141. Here 
we check the validity of the |Paardekooper et al.| 
(2011) formula in these more realistic discs. 

We consider the case of an accretion disc 
with constant M = 410 _8 MQ/yr and we com¬ 
pute the total torque for planets in the mass in¬ 
terval 20 — 70 Earth masses placed at distances 
of 4AU, 5AU and 6AU from the central star. We 
have also made a few tests on smaller planetary 
masses. The values of the masses and resolutions 
(N r , Ng, N v ) for the SAD disc are shown in Ta¬ 
ble [2] 


We have used a resolution of 7 grid-cells in 
the half width of the horseshoe region for all 
the computations in the interval 20 — 70 Earth 
masses. We have checked that this is needed in 
order to the have results independent of the res¬ 
olution. For small planetary masses a resolution 
of 4 grid-cells is enough to achieve convergence 
(see table |2|). 

In Fig|5| (top) we report the total torque 
computed on planets held on fixed orbits at 
4AU. The transition from inwards to outwards 
migration of small masses occurs for planet 
masses about twice as big as predicted by the 
Paardekooper et al. formula. We find again the 
“cold finger” effects, with a transition to negative 
torque at about 10 Earth masses. 

In the interval [15,3O]M0 we observe a 
positive total torque with values smaller than 
expected from the formula as in the previ¬ 
ously examined RED and SED. In the inter¬ 
val [40, 7O]M0 the results from numerical sim¬ 
ulations very nicely agree with the torque pro¬ 
vided by the analytical formula. We find positive 
torques up to planetary masses of 6OM0. 

Similar results are found at 5AU and 6AU 
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mass (Mq) (N r ,Ng, N v ) n cells in Xf, a 2 n distance (AU) 


5 

(1600,110,1132) 

4 

4 

15 

(1320,96,910) 

4 

4 

15 

(1256,90,886) 

6 

5 

20 

(1066,74,750) 

4 

4 

20 

(1600,100,1132) 

7 

4 

20 

(1318,80,920) 

6 

5 

20 

(1528,90,1072) 

7 

5 

20 

(1134,90,758) 

7 

6 

30 

(1502,90,1068) 

7 

4 

30 

(1250,76,878) 

7 

5 

30 

(1080,78,718) 

7 

6 

40 

(1360,86,972) 

7 

4 

40 

(1080,70,758) 

7 

5 

40 

(926,66,622) 

7 

6 

50 

(1226,76,868) 

7 

4 

50 

(966, 70, 680) 

7 

5 

50 

(826,66,558) 

7 

6 

60 

(1118,80,794) 

7 

4 

60 

(882, 66, 624) 

7 

5 

60 

(754,66,514) 

7 

6 

70 

(1034,74,730) 

7 

4 

70 

(808, 66, 576) 

7 

5 

70 

(692,56,476) 

7 

6 


Table 2. Simulations parameters for set SAD. 


(Fig[5j middle and bottom panels), precisely, 
cores of 5OM0 undergo outwards migration at 
5AU and at 6AU. From the migration map the 
transition between outwards and inwards migra¬ 
tion occurs for slightly larger masses; however 
the overall picture and the radial extent of the 
outwards migration is quite well reproduced. 

As we pointed out in Section 3 this is an 
important result for models of giant planet for¬ 
mation. This may solve the problem pointed out 


Coleman & Nelson (2014) where all planets 


were lost by migration before becoming gas-giant 
planets. However, we stress that this is true only 
for discs with large M like the one studied here. 
As pointed out in Bitsch et al. (2014 2015 I, the 
maximum planet mass for outwards migration 
decreases with decreasing M. 


5.1 Radial torque distribution 

To further confirm that the outwards migration 
detected in the accretion disc corresponds to that 
previously found in radiative discs without stel¬ 
lar heating and without a net mass flux through 
the disc we show the contribution to the torque 
in the close vicinity of the planets for both SAD 
and RED. We compute the radial torque density 
F(r), i.e the torque exerted on the planet by a 
ring of disk material located at a distance r from 


the star. The integral of F(r) on the radial co¬ 
ordinate provides the total torque F to t shown in 
the previous section: 


Ttot — 


nr 
J r„ 


r(r)dr . 


( 8 ) 


Fig[6] shows the radial distribution of the torque 
exerted on planets of different masses in the ac¬ 
cretion disc with M = 4 10 -8 A Iq/x/v at 5AU (top 
panel) and the radial torque density for the RED 
disc for planets placed at 5.2AU (bottom panel). 

In a disc without thermal effects (isother¬ 
mal disc) r(r) would be positive for r < r p and 
negative for r > r p . Moreover, F(r)/Fo would 
be independent of the planet’s mass. Here we 
observe that r(r)/Fo changes with the plane¬ 
tary mass. This effect is due to the saturation of 
the entropy-driven corotation torque ( Bandeau] 
fe Masset|2008 Kley et al.|2009 1 which is a func¬ 
tion of planetary masses. We remark in Fig|6] 
that, when the planetary mass increases from 20 
to 70M®, the absolute value of F/ro decreases 
and the location where T transits from positive 
to negative, which is at r > r p for the 20 
planet, approaches r p . The displacement of the 
point where T(r) = 0 relative to r p is indica¬ 
tive of the strength of the total positive torque. 
In this case, the total torque decreases for in¬ 
creasing planetary mass, until it becomes slightly 
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0 10 20 30 40 50 60 70 

M/MEarth 



M/MEarth 

Figure 5. Accretion disc: comparison with 
Paardekooper et al. (2011) formula. Top panel: 
planets at 4AU, middle: planets at 5AU bottom 
panel: planets at 6AU. The red points are obtained 
with a resolution of 7 grid cells in the half width 
of the horseshoe region, the green ones with a 
resolution of 4 grid cells. 




Figure 6. Radial torque density for planet of dif¬ 
ferent masses at 5AU in the accretion disc with 
M = 410 ~ s M@/yr (top) and in the RED disc for 
planets at 5.2AU (bottom). 

negative for the 70M© planet. In the case of the 
RED disc the results are very similar, although 
one can notice quantitative differences. In this 
case, E(r p ) = 0 for a planet of 40Af®. In fact, the 
transition to inwards migration occurs at about 
40M®, as shown in Fig. [4] 


6 ON THE WIDTH OF THE 
HORSESHOE REGION 

We now go beyond the raw comparison of the 
torques predicted in Paardekooper et al. formula 
and measured numerically, in order to under¬ 
stand how such a good agreement is achieved 
for large masses but not for low masses. This 
is indeed surprising because the formula is de¬ 
rived in principle in the linear regime, while for 
m p > 3OM0, where the agreement is best, the 
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non-linearities should start to appear (because of 
the condition q > h 3 ). Here we focus on the size 
of the corotation zone, which is a key parameter 
for the estimate of the corotation torque and of 
the torque saturation. 

According to Masset et al. ( 2006| ) Eq|3j ob¬ 
tained for 2D discs, is valid only for planet to 
star mass ratio q < h 3 . The horseshoe region 


for large masses (q > 1.5 10 4 in Masset et al. 


(20061) behaves as in the restricted three body 
problem (RTBP, hereafter), i.e. scales with q 1 ' 3 . 
For intermediate masses, the width of the horse¬ 
shoe region is larger than what predicted by the 
q 1 ' 2 scaling. This is a manifestation of the flow 
nonlinearity turning out in a boost of the coro¬ 
tation torque (Masset et al. (20061). 

In this section, we provide a measure of the 
half width of the horseshoe region, ( Xhs 3 D here¬ 
after), for SED and for the SAD with planets 
respectively at 4AU and at 5AU and compare it 
to Xhs2D of Eq0 Since differences with respect 
to Eqj3]can come from 3D effects as well as from 
radiative effects we will proceed in 2 steps: i) we 
determine the difference between the horseshoe 
width in 2D and 3D simulations in an isothermal 
setting, ii) we compare isothermal 3D to fully ra¬ 
diative 3D. 


0.044, this value is used in the plot of the Xh S 2 D 
law. 

We remark that in 2D the effect of the soft¬ 
ening is felt well outside the Hill radius, so that 
the horseshoe region clearly depends on the soft¬ 
ening length. In all our 3D numerical simulations 
we have used a softening length of 0.6 Rh- In or¬ 
der to compare 2D and 3D isothermal results we 
have also run a set of 2D simulations with a soft¬ 
ening length of e = 0.6 H. In FigjT] we observe 
that in this case the transition between the q 1/12 
and the q 1 '' 3 scaling occurs in the mass inter¬ 
val 910~ 5 < q < 1.65 1CT 4 . The width of the 
horseshoe region, both in the linear and in the 
RTBP regime, is reduced of about 10% with re¬ 
spect to the values obtained for the smaller soft¬ 
ening length. In Fig[7]we have also reported the 
measure of the width of the horseshoe region 
obtained with 3D isothermal simulations with 
e = 0.6 Rh- The disc parameters are the ones 
of SAD set, with planets at 5AU. The aspectra- 
tio at this distance from the star is h = 0.044 
(FigJTJ). For 3D simulations the measure of the 
width of the horseshoe region is determined by 
computing the streamlines on the disc midplane. 
Results are in good agreement with the 2D case 
with e = 0.6 H 0 so that we do not observe 2D 
versus 3D effects in the horseshoe region width. 


6.1 Comparison between 2D and 3D 
simulations 


For the two dimensional case we used the 


FARGO code (Masset (2000)). We recall that in 
2D a softening length e is applied to the planet 
potential through: 


$ = - 


Gm v 


\Ja% + e 2 


(9) 


In 2 dimensional models, Eq(9] allows to mimic 
the average influence that the planet would have 
on the vertical gas column. The measure of 
the width of the horseshoe region is determined 
by computing the streamlines. In Fig[7]we plot 
the half width of the horseshoe region normal¬ 
ized over the planet semi-major axis a p for dif¬ 
ferent values of the planet to star mass ratio 
q, in the mass range considered in this paper 
([10, 70 Using the nominal setting given in 
|Masset et al.| (20061 we recover their results (see 
their Fig.9). In Fig 7]the half width of the horse¬ 
shoe region obtained with a softening length of 


e = 0.3 H (as in Masset et al. (20061) is fitted by 


xM h s2D — 2.45a(g/3) 1/3 


( 10 ) 


for q > 1.35 10 4 while the Xhs 2 D law of Eq| 
nicely fits the data for q ^ 4.5 x 10 5 . The value 
of the aspectratio for these simulations is h = 


6.2 Comparison between 3D isothermal 
and 3D radiative simulations 

We call xMhs 3 D — 0 .9xMhs2D the width of 
the horseshoe region obtained by the fit of 3D 
isothermal simulations in the RTBP regime and 
Xhs 3 D — 0.9xhs2D the fit obtained for the linear 
regime. Fig[8] shows the measure of the width of 
horseshoe region for SED with planets at 4AU0 
and SAD with planets at 5AU. 

For comparison we plot the results of the 
3D isothermal simulations of Fig[7] We recover 
a regime with a q 1 ^ 2 scaling and a regime with 
a q 1 ^ 3 scaling and in both regimes the horse¬ 
shoe region is narrower than in the isothermal 
case of about a factor 7~ 4 / 4 as expected from 


Paardekooper & Papaloizou (2009). In the same 
Fig[8]we have plotted the law for the horseshoe 
region used in the Paardekooper formula, i.e. 

x h s2Dl~ 1/A - 


5 A value of the softening length of e = 0.7 H, is 
shown to provide a radial torque density very similar 
to the one obtained with 3D simulations (seejKley et| 

|aJ7] ([2012]).) 

° The aspectratio for SED at 4AU is h = 0.041 
(FigQ. To be compared to the other data sets having 


h = 0.044 a rescaling by \ h/h is applied. 
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Figure 7. Simulations with isothermal EOS. Half¬ 
width of the horseshoe region normalized over a p as 
a function of the planet to star mass ratio q. Results 
are shown for two sets of 2D simulations with respec¬ 
tively e = 0.3 H and e = 0.6H and for 3D simulation 
with isothermal equation of state and e = 0 .6Rh- 
The lines show the different fits in q 1 ' 2 ( with labels 
%hs 2 D and 0-9xh S 2D, see text) and g 1 / 3 (with labels 
xMhs 2 D an d 0.9xMh S 2D » see text). 



Figure 8. Same as Fig[7] for 3D simulations for 
respectively SED with planets at 4AU, and SAD 
with planets at 5AU and for the 3D isothermal 
run of Fig@ The lines show the fits in q 1 / 2 and 
q 1 / 3 . The law for the horseshoe region used in the 
Paardekooper formula, 1 ' 4 , is also shown. 


We observe that the measured horseshoe 
width is respectively ~ 10 % narrower than 

Xhs 2 D , y~ 1 ^ A in the linear (or q 1 ^ 2 scaling) regime 
and about ~ 10 % larger in the nonlinear (or q 3 
scaling) regime. Now a larger horseshoe region 


causes saturation earlier. We discuss this point 
in Section 6.4. 

In the isothermal case the transition be¬ 
tween the two regimes is associated to a boost 
in the corotation torque ( Masset et al.|2 006) and 
a small softening length is required to detect it. 
Therefore, we can wonder if we have missed some 
effects in the torque computation by using a soft¬ 
ening length of 0.6Rh- 


6.3 Comparison between 3D radiative 
discs with different softening length 

In the case of 3D isothermal simulations sim¬ 
ulations it has been shown ( |Kley et al.| [2009) 
that a decrease in the softening length for a 
planet of 20 corresponds to a drastic change 
in the density structure in the planet vicinity 
which has an important impact on the corota¬ 
tion torque. Precisely, the authors found a torque 
excess in agreement with | Masset et al.| ( |2006[ |. 
The situation changes when taking into account 
thermal effects. In this case a deeper poten¬ 
tial corresponds to an increase of temperature 
near the planet so that the density distribution 
close to the planet is not drastically changed. 
In jKley et ah . (2009) (their Fig.14, top ) chang¬ 
ing the smoothing length from e = 0.8Rh and 
e = 0.5 Rh in a fully radiative setting acts in 
increasing the torque of about 5%, with respect 
to 30% of increase in the corresponding isother¬ 
mal case. We have extended the discussion of 
(Kle y et al.|2009| ) concerning the dependence of 
the torque on the softening length by comput¬ 
ing the width of the horseshoe region and the 
total torque acting on fixed planets in the range 
[10 : 6 O]M 0 for the RED set with t = 0.3 Rh- 
In Fig[9] it appears clearly that the width of the 
horseshoe region is almost not affected by the 
change in the softening length]^] The measured 
torque is also almost not affected by the change 
in the softening length. 

We observe that in the q 1 ^ 2 regime the 
horseshoe region is nicely fitted by the law used 
in the Paardekooper formula (Fi g§- However, 
although the measured width of the horseshoe 
region perfectly matches the law used in the for¬ 
mula up to 2 OM 0 the measured torque fits only 
qualitatively the one provided by the formula as 
explained in Section 4.2. This is not surprising 
since many other parameters enter in the for¬ 
mula and can possibly be different in realistic 
3D discs. 


' The aspectratio for RED at 5.2AU is h = 0.04 
(Fig[T]). To be compared to data having h = 0.044 a 


rescaling by \ h/h is applied. 
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Figure 9. Half-width of the horseshoe region as a 
function planet to star mass ratio ( q ) for 3D sim¬ 
ulations of RED set, for respectively e = 0.5 h'// 
and e = 0.3 Rh- The lines show the fits in q 1 / 2 
and g 1 / 3 obtained for the stellar irradiated discs 
(Fig[8]|. The law for the horseshoe region used in the 
Paardekooper formula, Xf la 2 n'Y~ 1//4 , is a l so shown. 


In isothermal discs a torque excess is cor¬ 
related to the departure from the linear regime 
(Masset et al. 20061 while in all our radiative 
discs we do not observe the same phenomenon. 
Taking into account thermal effects, the tran¬ 
sition from the linear to the nonlinear regime, 
even considering a small softening length (Fig[9|, 
is smoother in radiative discs (Fig[8] and Figj9| 
with respect the isothermal case with small soft¬ 
ening length (Fig0. Moreover, the mass at 
which we observe a departure from the linear 
regime in Fig[8]and Fig[9]corresponds to the on¬ 
set of torque saturation in Fig® and [5] middle 
panel. 


6.4 Impact of the measured horseshoe 
width on planet migration 

The width of the horseshoe region is of crucial 
importance for the saturation of the torque act¬ 
ing on more massive planets. A larger width 
of the horseshoe region increases the viscous 
timescale (Eq. [6j making torque saturation eas¬ 
ier. As demonstrated in Fig[8] and Figj9j the 
width of the horseshoe region increases for more 
massive planets and does not follow the simple 
law of Xhs 2 D , y —1 , which is only valid for low 

mass planets. 

I 11 Fig |10| we show the migration map of 
the SAD disc, where we used the width of the 
horseshoe region measured by our simulations 
(Fig|5| instead of X/ lS 2 o 7 -1 ^ 4 - When compar¬ 


Figure 10. Migration maps for simulation set SAD 
with M = 4 10“ TIq/ yr obtained using the width 
of the horseshoe region measured by our simulations 
(Fig[sJ) instead of Xhs 2 Dl 1 / 4 . 


ing to Fig[2jtop panel we see that for plane¬ 
tary masses of up to 4OAI0 the migration map 
remains unchanged while for larger planetary 
masses, we observe that the transition of the 
torque from positive to negative occurs in the 
interval 40 — 5OM0. The transition occurs in the 
interval 60 — 7OA/0 in Fig[2]top panel. This is 
caused by a larger width of the horseshoe re¬ 
gion, which makes saturation of the corotation 
torque easier. Hence the total torque transitions 
into negative values at smaller planetary masses. 
Adapting the width of the horseshoe region to 
our measured values in the torque formula there¬ 
fore over predicts the transition to inward migra¬ 
tion compared to our simulations. 

The difference between our simulations and 
the torque predicted by the formula is therefore 
not only caused by a different width of the horse¬ 
shoe region. The cooling process plays also an 
important role and the vertical cooling in 2D disc 
is treated like blac.kbody radiation, while our 3D 
disc features vertical heat diffusion. 

Therefore we think that the quantita¬ 
tive agreement between the torque formula of 


Paardekooper et al. (2011) and simulations for 


planetary masses mp > 30 Mq is a coincidence 
and the formula should be updated with results 
of 3D simulations. However, this is beyond the 
scope of this paper and left for future work. 


7 CONCLUSION 

Using 3D hydrodynamical simulations including 
stellar and viscous heating as well as radiative 
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cooling we have computed the torque acting on 
planets of various masses kept on fixed circular 
orbits in both equilibrium and accretional pro¬ 
toplanetary disks. We confirm the results previ¬ 
ously obtained on constant viscosity equilibrium 
discs, where the only source of heating was pro¬ 
vided by viscous friction. The comparison be¬ 
tween the total torque obtained with numeri¬ 
cal simulations and the one predicted by semi- 
empirical formula (Paardekooper et al. 2011 ) 
is in quantitatively good agreement for large 
masses (> 2 OM 0 ). In particular, the formula 
predicts well the maximum planet mass at which 
inwards migration is prevented by the action of 
the entropy-driven corotation torque. 

Instead, although our results remain similar 
to those expected by the Paardekooper formula, 
we observe quantitative differences as large as 
a factor of ~ 2 concerning: (i) the value of the 
maximum positive torque as a function of the 
planet’s mass, (ii) the value of the planet mass 
at which the torque is maximal and (iii) the min¬ 
imum planet mass allowing outwards migration. 
The last difference is due to the presence of a 
’’cold finger effect” already discussed extensively 


Lega et al. (20141. The differences concerning 


the items (i) and (ii) have already been pointed 
out by Kley et al. (2009) as caused by differ¬ 
ences between 2D and 3D effects. Moreover, by 
measuring the width of the horseshoe region we 
have observed that the value of the planet mass 
at which the torque is maximal is correlated to 
the departure from the linear regime. For large 
planetary masses the quantitatively good agree¬ 
ment between the numerical results and the pre¬ 
diction of the analytic formula was unexpected 
since the formula is based on unperturbed discs 
while planets with masses larger than 2 OM 0 
start opening a partial gap. Precisely, torque 
saturation has the same behavior in the torque 
formula and in simulations though saturations 
effects appear for larger planetary masses with 
respect to the torque formula. From the mea¬ 
sure of the width of the horseshoe region, we 
have found that the onset of saturation is cor¬ 
related to the onset of the nonlinear regime not 
taken into account by the formula. Plugging the 
measured horseshoe width in the formula gives a 
worse comparison. Therefore, we think that the 
quantitative agreement on large masses is a co¬ 
incidence, and future works is needed to provide 
a more accurate torque expression. 

At the present state of the art, though there 
are some differences between the torque formula 
and 3D simulations, the zero torque location at 
the transition from outwards to inwards migra¬ 
tion at large planetary masses is reproduced sur¬ 
prisingly well, making the formula applicable to 


the study of planetary migration properties from 
the unperturbed disc structure. 
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